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Abstract 

The objective of this paper is to study the dynamical behaviour systematically of an 
ecological system with Beddington-DeAngelis functional response which avoids the criticism 
occurred in the case of ratio-dependent functional response at the low population density 
of both the species. The essential mathematical features of the present model have been 
analyzed thoroughly in terms of the local and the global stability and the bifurcations arising 
in some selected situations as well. The threshold values for some parameters indicating the 
feasibility and the stability conditions of some equilibria are also determined. We show 
that the dynamics outcome of the interaction among the species are much sensitive to the 
system parameters and initial population volume. The ranges of the significant parameters 
under which the system admits a Hopf bifurcation are investigated. The explicit formulae 
for determining the stability, direction and other properties of bifurcating periodic solutions 
are also derived with the use of both the normal form and the central manifold theory (cf. 

Carr [1]). Numerical illustrations are performed finally in order to validate the applicability 
of the model under consideration. 

Mathematics Subject Classification: 92D25, 92D30, 92D40. 
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1 Introduction 

Mathematical model is an important tool in analyzing the ecological models. Ecological prob¬ 
lems are challenging and important issues from both the ecological and the mathematical point 
of view (cf. Anderson and May [2] , Beretta and Kuang [3] , Freedman [4] , Hadeler and Freedman 
[5], Hethcote et al. [6], Ma and Takeuchi [7], Venturino [8], Xiao and Chen [9]). The dynamic 
relationship between predator and its prey has long been and will continue to be one of the 
dominant themes in both ecology and mathematical ecology due to its universal existence and 
importance. The most common method of modelling that ecological interactions consists of 
two differential equations with simple correspondence between the consumption of prey by the 
admissible predator and their population growth. The traditional predator-prey models have 
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been studied extensively (cf. Cantrell and Cosner [10], Cosner et al. [11], Cui and Takeuchi 
[12], Huo et al. [13] and Hwang [14]), but those are questioned by several biologists. The 
most crucial element in these models is the “functional response” - the expression that describes 
the rate at which the number of prey consumed by a predator. Modifications were limited to 
replacing the Malthusian growth function, the predator per capita consumption of prey func¬ 
tions such as Holling type I, II, III functional responses or density dependent mortality rates. 
These functional responses depend only on the prey volume x, but soon it became clear that 
the predator volume y can influence this function by direct interference while searching or by 
pseudo interference (cf. Curds and Cockburn [15], Hassell and Varley [16] and Salt [17]). A 
simple way of incorporating predator dependence in the functional response was proposed by 
Arditi and Ginzburg [18], who considered this response function as a function of the ratio x/y. 
The ratio-dependent response function produces richer dynamics than that of all the Holling 
types responses, but it is often criticized that the paradox occurred at the low densities of both 
populations size. Normally one would expect that the population growth rate decrease when 
both the populations fall bellow some critical volume, because food-searching effort becomes 
very high. For some ecological interaction ratio-dependent model give the negative feed back. 
Thus, the Lotka-Volterra type predator-prey model with the Beddington-DeAngelis functional 
response has been proposed and well studied. Keeping these in mind, the proposed model can 
be expressed as follows: 


dx'^ 
dt 
dx'o 

~dr 

with the initial conditions x'^(O) = > 0 and X 2 ( 0 ) = x'.^ > 0. The functions x'i{t), x^it) are 

the volumes of prey and predator at any time t. All the system parameters are assumed to be 
positive and have their usual biological meanings. The functional response in system 

(1.1) was introduced by Beddington [19] and DeAngelis et al. [20] as a solution of the observed 
problem in the classic predator-prey theory. It is similar to the well-known Holling type-H 
functional response but has an extra term 61 X 2 in the denominator which models mutual inter¬ 
ference between predators. It represents the most qualitative features of the ratio-dependent 
models, but avoids the “low-densities problem”, which usually the source of controversy. It 
can be derived mechanistically from considerations of time utilization (cf. Beddington [19]) or 
spatial limits on predation. 

The present study under consideration has been carried out sequentially in the latter sections 
as follows: The basic assumptions and the model formation are proposed in Section 2. Section 
3 deals with some preliminary results. The equilibria and their feasibility are rightly given in 
Section 4. The local analyses of the system around the boundary as well as interior equilibria 
are discussed in Section 5. The global analysis of the system around the interior equilibrium 
is studied at length in Section 6 . Simulation results are reported in Section 7 while a final 
discussion and interpretation of the results of the present study in ecological terms are rightly 
included in the concluding Section 8 . 




C\X-^X2 


ai-\-x'-^+hiX2 
cieix'-^X2 


- -61X2 + 


( 1 . 1 ) 
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2 Model formulation 


Firstly we replaced the logistics growth function rxi(l — of the prey species by the modified 
quasi-linear growth function rxi{l — ~ ^ order to make the 

model free from any axial equilibrium. Which fits better for some special type of ecosystem, 
whereof environmental carrying capacity varies w.r.t. its prey volume, i.e., carrying capacity 
is always greater than its present prey volume. In the present model we introduce one more 
predator species in the model (1.1) to make it one step closure to reality. Thus, our final model 
is extended to the following form: 


+ 

^ + 




CIXIX 2 


xi-\-k' ai-\-xi-\-biX2 
cieixiX2 

ai+xi-\-biX2 ’ 

0262x1x3 

a2+xi-\-b2X3 


C 2 X 1 X 3 

a2-\-xi-\-b2X3 


( 2 . 1 ) 


where xi is the population volume of the two prey species and X 2 , X 3 are the population 
volumes of the predator species at any time t. It is assumed that all the system parameters 
are positive constants. Here r and k are the growth rate and the half-saturation constant for 
the prey species, di, <52 are the first and second predators death rate respectively, ci, C 2 are 
the respective search rates of the first and second predator on the prey species, ^ are the 
maximum number of prey that can be eaten by the first and second predator per unit time 
respectively; ^ being their respective half saturation rates while ei, 62 are the conversion 
factors, denoting the number of newly born first and second predator for each captured prey 
species respectively (0 < ei, 62 < 1). The parameters bi and 62 measure the coefficients 
of mutual interference among the first and second predator species respectively. The terms 
ai+^xi+bix 2 a 2 +Ti^+b 2 X 3 denote the respective predator responses on the hrst and second prey 
species. This type of predator response function is known as Beddington-DeAngelis response 
function (cf. Beddington [19] and DeAngelis et al. [20]). 


3 Some preliminary results 

3.1 Existence and positive invariance 

Letting, x = (xi,X 2 ,X 3 )^ / : —)• R^, F = {fi, f 2 , fsY, the system (2.1) can be rewritten as 

X = f{x). Here fi G C°°(R) for i = 1,2,3, where /i = rxi{l-^) - a,lZ+Yx 2 " 

/2 = -hx 2 + aYi 7 Y+biX 2 -^3 = -hx 3 + ^Y+YY+b^xs ' S^ce the vector function / is a smooth 
function of the variables {xi,X 2 ,X 3 ) in the positive octant 0 ° = {(xi,X 2 ,X 3 ) : xi > 0 ,X 2 > 
0 , X 3 > 0 }, the local existence and uniqueness of the solution of the system ( 2 . 1 ) hold. 

3.2 Persistence 

If a compact set D C = {(xi,X 2 ,X 3 ) : Xi > 0, i = 1,2,3} exists such that all solutions of 
(2.1) eventually enter and remain in D, the system is called persistent. 

Proposition 3.1. The system (2.1) is persistent if the conditions: {i)r > (5i + 82 , (n) xij > 
{iii)xi 2 > are satisfied. 
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Proof. We use the method of average Lyapunov function (cf. Card and Halam [21]), considering 
a function of the form 

V{xi,X2,X3) = xj^x'l^x'l^, 

where 71,72 and 73 are positive constants to be determined. We define 


n(xi,X 2 ,X 3 ) 


V 


71 r- 


C 1 X 2 


+73 (-^2 + 


rxi 

xi + k ai + xi + 61 X 2 
€ 262 X 1 


a2 + xi + 62 X 3 


£ 2^:3 \ ^ J 5 ^ \ cieixi 

02 + 3^1 + 62X3/ ^ V ai+xi + 6iX2 


We now prove that this function is positive at each boundary equilibrium. Let 7 ^ = 7 , for 
i = 1, 2, 3. In fact at Eq, we have 11(0, 0,0) = r — 61 — 62 > 0 from the condition (i). Moreover, 
from condition (n) and {Hi), we find the values of 11 at Ei and E 2 respectively. 


n(xi„X2„0)=7(-52 + ^^5^) >0, 
n(xi2,o,x32) = 7(-(5 i + > 0. 

V oi + X 12 / 


Hence, there always exists a positive number 7 such that H > 0 at the boundary equilibria. 
Hence V is an average Lyapunov function and thus, the system (2.1) is persistent. □ 

Since the system is uniformly persistent, there exists cr > 0 and r > 0 such that Xj(t) > cr, for 
all t > T, i = 1, 2, 3. 


3.3 Boundedness 


Boundedness implies that the system is consistent with biological significance. The following 
propositions ensure the boundedness of the system ( 2 . 1 ). 

Proposition 3.2. The prey population is always bounded from above. 


Proof. Before proving that the prey population is bounded above, we need to prove that the 
predator populations X 2 and X 3 are bounded above. To prove this result, considering the second 
sub equation of the system ( 2 . 1 ) and one can obtain the following differential inequality: 


-w - 


ciei)x 2 . 


Integrating the above differential inequality between the limits 0 and f, we have X 2 {t) < 
X 2 ( 0 )e“^^i“‘^i®i)h Thus, if (61 — cici) > 0, then it is obviously found a positive number ri 
there exists a positive constant mi such that X 2 {t) < mi, for all t > ti. By using the similar 
argument, one can obtain that, if {62 — 0262 ) > 0 , then corresponding to a positive number T 2 
there exists a positive constant m 2 such that X 3 (t) < m 2 , for all t > T 2 . Both the results can 
be written unitedly as Xi > m = min (mi, m 2 ) for all t > T 3 = max(ri,r 2 ), i = 2,3, with the 
additional condition min (61 — cici, 62 — 0262 ) > 0 . 
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Now from the first sub-equation of ( 2 . 1 ), the following inequality is found 


dxi ^ ((ei -h 62)0- - rk)xi /k[rm - (ei -h 62)0-) 
dt ~ km V (ei -|- 62)0- — rk 


Xl 


Hence, by using standard but simple argument, we have 


krm — (ei -|- e2)ka rm rk 

limsupa:i(t) < —-r-— = w, where - < a < -. 

t^+00 [ei + e2)a-rk ei + 62 ei + 62 


□ 

Proposition 3.3. The solutions of ( 2 . 1 ) starting in are uniformly bounded with an ultimate 
bound. 

Proof. Considering the total environment population x = xi + ^ + ^. Using the theorem 
on differential inequality (cf. Birkhoff and Rota [ 22 ]) and following the steps of Haque and 
Venturino [ 23 ], Sarwardi et al [ 24 ], boundedness of the solution trajectories of this model is 
established. In particular, 

limsup (xi + — + —)< (^ + + = M, where p = min ( 1 ,( 5 i, 52 ), ( 3 . 1 ) 

t^+00 ^ Cl 62 / p 

with the last bound is independent of the initial condition. 

Hence, all the solutions of ( 2 . 1 ) starting in for any 6 > 0 evolve with respect to time in the 
compact region 

H = I {xi,X2,X3) e Ri — -I-— ( 3 . 2 ) 
I ei 62 J 


4 Equilibria and their feasibility 


The equilibria of the dynamical system ( 2 . 1 ) are given as follows: 

1 . The trivial equilibrium point i?o( 0 , 0 , 0 ) is always feasible. 


2 . (a) The first boundary equilibrium point is £'i(xii,X 2 i, 0 ). The component is a root 
of the quadratic equation hx^^ + {I2 + hk + rfc 6 i 6 i)xi^ -|- 12k = 0 , where li = (cii — ci 6 i), 
I2 = a\ 5 i. If h < 0 , then the quadratic equation in xij possesses a unique positive root and 


consequently X2i 


(ciei-(5i)a:ij^-ai(5i 

61 5 i 


. The feasibility of the equilibrium Ei is maintained if the 


condition xi, > is satisfied. 

-*^1 c\e\—o\ 


2. (b) The second boundary equilibrium point is £'2(2:12, 0 ,X32). The component X12 is the root 

of the quadratic equation mix\^ + {m2 + mik + rkh2e2)xi.^ + m2k = 0, where mi = (^2 — C262), 

I2 = 0252 - If "2-1 < 0 , then the quadratic equation in xij possesses a unique positive root and 

consequently xsj = *^62^2^^ —^2^ ^ feasibility of the equilibrium £2 is maintained if the 

condition xi, > holds. 

^2 ^ 262—02 


3 . The interior equilibrium point is £*(xi*,X2*,X3*), where the first component xi* is the root 
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of the following quadratic equation: 


+ (n 2 + nik + rkbib2eie2)xu + n2k = 0 , ( 4 . 1 ) 

where ni = biei{52 - 0262) + 6262(51 - cici) and 77-2 = 62 e 25 iai + 6 iei 52 a 2 - 

Case I: Let ni < 0 . In this case there exists exactly one positive root of the quadratic equation 
( 4 . 1 ) irrespective of the sign of (122 + nik + rkbib2eie2)- 

Case II: Let ni > 0 . In this case there are two possibilities: (i) if n2 + nik + rkbib2eie2 > 0 , 
then there is no positive solution and (ii) if n2 + nik + rkbib2eie2 < 0 , then there exists two 
positive roots or no positive root. 

In this present analysis we consider the Case 1 . Under this assumption the next two components 
of the interior equilibrium can be obtained as X2* = ^ = (C2e2-fe^)x^i*-a2^2 ^ 

The feasibility of this important equilibrium point is conhrmed under the condition xi* > 
max I I • Moreover, the positivity condition of second and third components of 

the interior equilibrium ensures the impossibility of the Case II. 

Remark: The feasibility and existences conditions of both the planer equilibria Ei and E2 
immediately implies the existence of the unique feasible interior equilibrium point . But the 
existence of the unique feasible interior equilibrium point E^ implies three possibilities: (i) Ei 
exists and E2 does not exist, (ii) E2 exists and Ei does not exist, (iii) existence of both. 


5 Local stability and bifurcation 

The Jacobian matrix J{x) of the system ( 2 . 1 ) at any point x = {xi,X2-,x^) is given by 

C2X\(a2+xi) 

(02+3:1+622:3)^ 

0 

r I £2622:1(02+3:1) 

2 ((U+Tl+fe+fF / 

( 5 . 1 ) 

tr(J), k2 = M and 

Note that the conditions for occurrence of Hopf bifurcation are that there exists a certain 
bifurcation parameter r = such that C'2(rc) = A:i(rc)/c2(u) ~ k^i^c) = 0 with k2{rc) > 0 and 
^(Re(A(r)))|r=rc / 0 , where A is root of the characteristic equation A(A) = 0 . 

5.1 Local analysis of the system around Eo, El, E 2 

Stability: The eigenvalues of the Jacobian matrix J{Eq) are r,—6i and —52- Hence Eq is 
unstable in nature (saddle point). Let J{Ei) = (Cij)3x3 and J(£^2) = (hii)3x3- Using the 
RouthHurwitz criterion, it can be easily shown that the eigenvalues of the matrices J{Ei) and 
J{E2) will have negative real parts iff the conditions eixi^ +X2^ > and 62X12 +X32 > 

k{i 6262) 02 j-ggpeci^ively. Hence the equilibria Ei and E2 are locally asymptotically stable under 


^(x) 3 x 3 = 


/ 

(Si+fcp’ 


V 


01x2(01+61x2) _ 62x3(02+62x3) 
(oi+X1+6iX2)^ (o2+X1+62X3)^ 

0101x2(01+61x2) 

(oi+xi+ 6 iX 2)^ 

6262X3(02+62X3) 

(o2+Xi+62X3)^ 


6ixi(oi+xi) 

(oi+xi+ 6 iX 2)^ 

_ e I 6ieixi(ai+xi) 
1 ' (oi+xi+6iX2)^ 

0 


Its characteristic equation is A(A) = A^ + /ciA^ + k2X + A:3 = 0 , where ki = - 
ks = — det( J); M being the sum of the principal minors of order two of J. 


6 



















the conditions eiXij +X2i > and 62X13 +X32 > respectively (cf. Section 

4 of of Sarwardi et al. [ 25 ]). 

Bifurcation: Since the equilibrium point Eq is a saddle in nature, hence, there is no question 
of Hopf bifurcation around this equilibrium. In order to have Hopf bifurcation around the 
equilibria Ei, E2, it is sufficient to show that the coefficient of A in the quadratic factor of the 
characteristic polynomial of J{Ek) {k = 1 , 2 ) is zero and the constant term is positive. The 
conditions for which annihilation of the linear terms in the quadratic factors of the characteristic 
polynomials of J{Ei) and J{E2) can be made possible are + ^22 = 0 and 7711 + 7733 = 0 . 
For detailed analysis, interested readers are referred to Appendix A of Haque and Venturion 
[ 26 ] . The parametric regions where Hopf bifurcations occur around Ei and E2 are respectively 
established by the equality constraints eixi^ +X2i = and 62X13 +X33 = ^262) 02 ^ 


5.2 Local analysis of the system around the interior equilibrium 


Proposition 5.1. The system ( 2 . 1 ) around E^ is loeally asymptotically stable if the eondition 
(i) k < min{ai + 61X2*, 02 + 62X3*} is satisfied. 


Proof. Let J(x*) = is the Jacobian matrix at the interior equilibrium point E^ = x* of 

, . ^ \ 7 cxxi,X2,{k-{ai+hiX2)) C2Xi,X3,{k-{a2+h2X3,)) 

the system ( 2 . 1 ). The components of J(x*) are Jn = (^,,^+k){a,+x,,+b,x2!) + {x,,+k)ia2+x,.,+b2xj) ^ 


cixit(ai+xit) 


C 2 Xit{a 2 +xit) 


_ cieiX 2 t,{ai+biX 2 t) 


•^12 - -(ai+xi*+&iX2*)^ < U, ^13 - - (a2+xi*+fe2a;3*)^ < *^21 “ (ai+xi*+6iX2.)" > "^22 “ 

bicieia;i*3;2* 7. _ _ ri T., — <^ 2623 : 3 *( 02 + 623 : 3 *) ^ q 7 n T — 62C2e23;i*3;3* 

“(ai+xi*+6i3:2*)2 < «^23 - U, J3I - (a2+3:i*+623:3*)2 > ‘J’ '^32 - U, J33 - “ (a2+3:i*+623:3*)2 < 

Then the characteristic equation of the Jacobian matrix J(x*) can be written as 


+ kiX^ + k 2 X k^ — 0, 


( 5 . 1 ) 


where ki = —tr(J) = —(Jn + J22 + J33), ^2 = Mu + M22 + M33 = (Jn J22 — J21J12) + J22J33 + 
(Jii J33 —J31 Jis): ks = —det (J) = —[J11J22J33 — J12J21J33 — J31J13J22), and C2 = kik2 — k3 = 
— (Jn + J 22 )(J 33 (Jll + J22 + J33) + (Jn J22 — J2IJ12)) + Jl 3 J 3 l(Jll + J 33 )- 

It is clear that /ci > 0 if Jn < 0 , i.e., k < min{ai + 61X2*, ^2 + ^2X3*} and consequently C2 > 0 . 
Hence the Routh-Hurwitz condition is satisfied for the matrix J*, i.e., all the characteristic roots 
of J* are with negative real parts. So the system is locally asymptotically stable around E^. 
Theorem 5 . 2 . The dynamical system ( 2 . 1 ) undergoes Hopf bifurcation around the interior 
equilibrium point E^ whenever the critical parameter value r = Vc contained in the domain 

Dhb = |xc G : C' 2 (rc) = ki{rc)k2irc) - ksivc) = 0 with A: 2 (xc) > 0 and o|. 


Proof. The equation (5.1) will have a pair of purely imaginary roots if kik2 — fcs = 0 for 
some set of values of the system parameters. Let us now suppose that r = Xc be the value of 
r satisfying the condition kik2 — k^ = 0 . Here only Jn contains r explicitly. So, we write the 
equation kik2 — k^ = 0 as an equation in Jn to find Xc as follows: 


h-i Jn + ^<2 Jn + /13 — 0 , 


( 5 . 2 ) 


where hi = J22 + J33, ^2 


- /2 
^22 


+ J? 


33 


Ji3 J31 — J12 J21! 


(J22 + J33) J22 J33 — Jl 3 J3IJ33 — 
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Jl2J2lJ22- 

Thus, Jii = ^(-^2 ± 4/11/13) = -^n- 

Or, 

_ (xu + r CiX2^,(ai + &lX2^,) 

(ai + xi* + 61X2*)^ 


02X3^ {a2 + fc 2 a^ 3 *) 

(02 + Xi* + b 2 X 3 ^Y 


= Tc- 


( 5 . 3 ) 


Using the condition kik2 — /cs = 0 , from equation ( 5 . 1 ) one can obtain 

(A + fci)(A^ + ^2) = 0 ) (^- 4 ) 

which has three roots Ai = +iy/k2, A2 = —iy/k2, A3 = —ki, so there is a pair of purely 
imaginary eigenvalues ±z\/fc2- For all values of A, the roots are, in general, of the form Ai(r) = 
?i(?’)+ *6(5"), A2(r) = 6 ( 5 ") - *6(*’), \3{r) =-ki{r). 

Differentiating the characteristic equation ( 5 . 1 ) w.r.t. r, we have 


dr 


Hence, 


y?k\ + Xk2 + ^3 I 
“ 3 A 2 + 2 A:iA + 1^2 ’ 
ks — 6^1 + 

2 {k 2 - ikl^/^) 

ks - {kik2 + kik2) .V^ikih + 66 - A:i66) 
2 (A;^ + 6) ^ 26 (A:i+ 6 ) 

^ + 666 feiV6^ - 

2(A:f + 6) I- 26 26(A:f + 6)- 


dr 


(Re(A(r)))|^^^^ 


dC 2 

_dr_ 

2 {kl + 6) 


r=rc 


^0. 


( 5 . 5 ) 


( 5 . 6 ) 


Using the monotonicity condition of the real part of the complex root *^(F 6 A(r))) ^ q 

Wiggins [ 27 ], pp. 380 ), one can easily establish the transversality condition ^ 0 , to 

ensure the existence of Hopf bifurcation around . 


6 Global analysis of the system around the interior equilibrium 

6.1 Direction of Hopf bifucation of the system (2.1) around 

In this Section we study on the direction of Hopf bifucation around the interior equilibrium. 
From the model equations ( 2 . 1 ), we have 


X = /(x). 


( 6 . 1 ) 
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where x = (xi,2:2,/ = (/^/^/^)* = 

Here, at x = x*, / = 0 . Let y = {yi,y2,yz) 
equation ( 6 . 1 ), we have 




ai+xi+hiX 2 

"1-^2 -r aij^x-i_+biX2 

SoX-i A _ '^ 2 e 2 XiX 3 

2 3^ a2A-x:iA-b2X3 

(Xi - Xi*, X 2 - X 2 *, X3 - X3*) 


£23:13:3 

a 2 A-xi+b 2 X 3 


Putting in 


y = J(X1*,X2*,X3*)?/ + 0, 


( 6 . 2 ) 


where the components of nonlinear vector function cj) = 02, 03)* are given by 


= fi 


,.2 , xi 


The coefficients of nonlinear terms in yi, i = 1 , 2 , 3 are given by 

= __2rfc2 2 eiX 2 (ai+fci 3 : 2 ) 2023:3(02+623:3) _ 26 icixi(ai+xi) ^1 _ 262023:1(02+3:1) 

■'3:13:1 (xi+fc)^ ' (01+3:1+61x2)^ ' (a2+xi+62X3)^ ’ 3X2X2 (ai+xi+6iX2)^ ’ •'X3X3 (a2+xi+62X3)^ ’ 

xl _ Q /•! _ oiai(oi+xi+ 6 iX 2 )+ 26 ioixiX 2 xl _ 02 O 2 (o 2 +xi+ 62 X 3 )+ 262 C 2 X 3 Xi , x 2 _ 

3 X 2 X 3 ~ 3 X 1 X 2 ~ (oi+Xi+ 6 iX 2 )^ ’ •' 3 : 33:1 ~ ( 02 +Xi+ 62 X 3 )^ ’ 3 x 1 X 1 “ 

20161x2(01+61x2) x2 _ 26 ioieixi(oi+xi) ^2 _ n f2 _ oiOiei(oi+xi+6iX2)+26i0161X1X2 

(ai+xi+ 6 iX 2 )^ ’ ■'3:2x2 ~ (oi+xi+ 6 iX 2 )^ ’ ■'3:33:3 “ ■'X1X2 ~ (oi+xi+ 6 iX 2 )^ ’ 

x2 _ n /■2 _ n. fS _ 20262x3(02+62x3) ^3 _ ri x 3 _ 2620262X1 (02+xi) 

3X2X3 — U, Jx^xi — U, Jxixi — (a2+xi+62X3)i* ’ ■'X2X2 ~ *3, jX3X3 “ (a2+xi+62X3)^ ’ 

x 3 =(] f 3 — O2 02 e 2 (o 2 + 3 i+ 62 X 3 )+ 262 0262x1x3 ^3 _ ^ 

3x2X3 *3 , ^3.33;^ (o2+Xi+62X3)^ ’ 3x1X2 

Let P be the matrix formed by the column vectors (u2,ui,U3), which are the eigenvec¬ 
tors corresponding to the eigenvalues Ai^2 = ±i\/^ and A3 = —ki of J(xii,,3:2*,a;3i,), then 
J(xi*,a:2*,a:3*)u2 = iV^u.2, J{xu,X2*,X3^)ui = and J(xi*,X2*,X3*)u3 = -kiu^. 


Thus, 


/ 0 1 

— J 2 lV^ —J21J22 

— —^ 31^33 

\ -^ 33+^2 '■^ 33+^2 


1 \ 

— J2I 

J 22 +kl 
-Jzi 
^ 33+^1 / 


{Pij)3x3- 


Let us make use of the transformation y = Pz, so as the system ( 6 . 2 ) is reduced to the following 
one 

/ 0 -iV^ 0 \ 

z = J(xi*,X2*,X3*)Pz-I-P “^0 = I 0 0 |z-|-P“^ 0 . ( 6 . 4 ) 

\ 0 0 -ki) 


Here P 


-1 _ Adjp 

~ detP 


(%)3x3, where 


911 

913 

922 

932 


_1_(_ J21J22J31 _ J21J33J31 _1 „ _ _1_ ( ,031 _ J31J33 \ 

detPl(j| 2 +fc 2 )(J 33 +A:i) (J|2+fc2)(J22+A:i)det P 1 J|3+fci 

1 ( -J2I I J21J22 „ _ 1 ( J 2 lVk^J 31 _ J 21 \/k 2 J 31 ] 

detpl J2i+fci I^^+^h 921 detpV(j|2+fc2)(J22+fci) (J|3+A:2)(J22+fci) 

1 \/fc2 J3I „ _ 1 —-y/faJai „ _ 1 / J21\/fc2 J3IJ33 _ J 2 lVk 2 J 3 lJ 22 \ 

detP(J|3+fc2)’ «23 detP^jf^T^: 931 det P 1 ’ 

1 -pk 2 J 21 „ _ 1 V^J 21 

det P (J|2+fc2) ’ det P ( 1 ^+^ ‘ 


The system ( 6 . 4 ) can be written as 



dZ 3 

dt 


-/ci 2;3 -kG(zi,2:2,2;3). 


( 6 . 5 ) 

( 6 . 6 ) 
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On the center-manifold (cf. Carr [ 1 ], Kar [ 28 ]) 


~ 2 + ^22^2) 


( 6 . 7 ) 


Therefore, 


0 -y/k^ fzi 


Z3 = ( +bl 2 Z 2 bi 2 Zi +b 22 Z 2 ) q j 

Using ( 6 . 4 ) and (6.6), we have 


— -\/^^122i + '\/^(^22 —&1i)2:i2;2 —■\/^^12^2 

( 6 . 8 ) 


is = -kiZ3 + q3i4>i + q32(t>2 + 93303 ( 6 . 9 ) 

From the equations (6.8) and ( 6 . 9 ), we have 

V^6i2^i + \/^(^22 — ^11)212:2 — y^bi 2 Z 2 

= -^kiibiizf + 2 bi 2 ZiZ 2 + b22zl) + 931 [fxixAPnZi + P12Z2 + Pi 3 -^{bnzl + 26122122 + 6222 !)}^ 
+^2222 +P 23 ^( 6 ii 2 i 26 i 22 i 22 + 62223)}^ + fAxalPSl^l + P 32 Z 2 + P 33 ^{bliz‘f 

+2bi2ZlZ2 + b22zl)A + 2/^33,JP31^1 +P32^;2 +P33^(6ii2i 26i22i22 + 62222 )}{pil 2 i + Pl222 

+Pi 3 ^( 6 ii 2 i + 26122122 + 62222)} + 2 fA^^{piizi +P12Z2 +Pi 3 ^( 6 ii 2 i -h 26122122 + 6222I)} 
'X{P 21 ZI +P22^2 +P 23 ^( 6 ll 2 i 26 i 22 i 22 + 6222I)}] 

+ 932 [/xixi{Pll 2 l +P12^;2 + 7 ' 13 ^( 6 ii 2 i 26 i 22 i 22 + 62223)}^ + + 7 ' 22^;2 

1 1 

+^232(^11^1 + 26122122 + 62222)}^ + +133222+^332(^11^1 + 26122122 + 6222!)}^ 


+2fl^^^{p2lZi +P22Z2 +P 23 ^( 6 ii 2 i 26122122 + 6222 |)}{p 3 i 2 i + P 32^2 +P 33 ^( 6 ll 2 i 


-F26i22i22 + 62222)} + 2/,;3,j,Jp3i2i + ^322:2 + P33 2 (^11^1 + 26l22l22 + 62222)}{pil2i + P12Z2 

+Pi 3 ^( 6 ii 2 i + 26122122 + 62222)} + 2/,^^3,2{pii2i + P12Z2 + Pi 3 ij{biizl + 26122122 + 6222I)} 

1 


x{P212:1 +P22^2 +P 237 j(^ll^l + 26 i 22 i 22 + 6222^)}] 


+933[/xia;i{l3ll2l +111222 + 1113^(^11^1 + 26l22l22 + 6222^)}^ + P22Z2 + P23 

X^(6 ii 2 i + 26122122 + 62222)}^ + +11322:2 +P33^(6ii2i 26l22l22 + 63223)}^ 

+ ‘^fl 2 X 3 {P‘il^l +P22Z2 +1123^(61121 + 26122122 + 6222 f)}{p 3 i 2 i +P 3222 +P 33 ^( 6 ii 2 i 
- h 26 i 22 i 22 + 62222)} + 2 /,^ 3 ,j,Jp 3 l 2 i +P 3222 +1133^(61121 + 26 i 22 i 22 + 62222)}{pil2i 
+P 12 Z 2 + Pi3]j{biizl + 26122122 + 6222I)} + 2fAx2{Pii^i +P 12 Z 2 +Pi3^{biizj 
+2bi2ZlZ2 + 6222|)}{p2l2l + P22Z2 + P23^{bnzf + 26i22i22 + 6222^)}] . 
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Comparing the coefficients of zf, Z\Z2 and z\ from both sides, we have 

'/^bi2 + ~^bu 

Q 3 I l^fxixiPll ffi fX 2 X 2 P 2 I ffi fxsx^Psi ffi fX2X3P‘21P^1 ffi “1“ 2/xia;2^^11^^2l] 

ffi932 \^fxixiPll ffi fX 2 X 2 P 2 I ffi fX 3 X 3 P 3 I ffi x 2X3P‘21P31 x3XiP3lPll “1“ 

+933 [fxixiPll + fx2X2P21 + fx3X3P31 + ^fx2X3P^lP3l + ‘^fx3XiP3lPn + ‘^fxiX2PnP2l] 

= Oi, (6.10) 


-'/h{b 11 — ^ 22 ) + ^ 1^12 

= ‘^Q3l[fxixiPnPl2 + fx 2 X 2 P 2 lP 22 + fx3X3P3lP32 + (P2lP32 + P22P31) + /igxi (PllP32 

+P 12 P 31 ) + /x\x2blll’22 +P 12 P 21 )] + 932[/xixiPllPl2 + fx2X2P2lP22 + /J 3 X 3 P 31 P 32 + fx2X3 
X{P21P32 + P 22 P 31 ) + /x3Xi(P11?’32 +P 12 P 31 ) + /xixa (^ 11^*22 +P 12 P 21 )] + 933 [/xixiPllPl 2 

+fx2X2P2lP22 + fx3X3P3lP32 + /xaxg (P2lP32 + P22P3l) + /xgxi (PllP32 + Pl2P3l) 
+/xiX2(PllP22 +PI 2 P 21 )] = ^ 2 , (6.11) 


and 


-2012 + ^^22 


-y^b 

931 [/xixil^l2 + fX 2 X 2 P 22 + fx3X3P32 + fX 2 X 3 P 22 P 32 + X 3 X 1 P 31 P 12 + 2/xiX2l^l2P22] 
+932 [/xixiPl2 + fx 2 X 2 P 22 + fx3X3P32 + ‘^fX 2 X 3 P 22 P 32 + fx3X\P3lPl2 + ^/xiX2Pl2P22] 

+933 [/xixil^l2 + fX 2 X 2 P 22 + fx3X3P32 + ‘^fX 2 X 3 P 22 P 32 + fx3X\P3lPl2 + 2/3;j^x2l^l2P22] 

= Os- (6.12) 


From equations ( 6 . 10 ), ( 6 . 11 ) and ( 6 . 12 ), we have 

\ki y% 0 

ki 

0 -\/fe ^ki 

The equation ( 6 . 13 ) gives the coefficients bu, 612 and 622 as follows: 

^2(f^i + f^3) ~ ^iV^^2 — kiQi) 



bu = 


bi 2 = 


b 22 = 


{^ + kik 2 ) 

^ _ fcl ^(03 - 111) 




y 4 + ^1^:2) 

fc 2 ( 0 l + 03 ) + ^ + ^^^2 
(|+A:ifc2) 


The flow of the central manifold is characterized by the reduced system as 


^ f Zl 

dt V 2^2 


0 -Vh 
0 


2:1 

22 


+ 


i?2 


( 6 . 13 ) 


( 6 . 14 ) 
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where = qu(i>i + qi 24>2 + qi 34>3 + h.o.t, = ^21^1 + q 224>2 + q 234’3 + h.o.t. The stability of 
the bifurcating limit cycle can be determined by the sign of the parametric expression 


n — Fj^ii + F112 + Fi2 


+ FL] 


where Fijk = origin. If the value of the above expression is negative, then the 

Hopf bifurcating limit cycle is stable and is called a supercritical Hopf bifurcation. If the value 
is positive, then the Hopf bifurcating limit cycle is unstable and the bifurcation is subcritical. 

Here 


-^11 xixiPll T fx2X2P21 T fx^x^Psi T fx^X]P3lP31 T fx\X2P^^P‘^^ ^912T fX 2 X 2 

XP2I + ‘^fx^X 2 PnP 2 l] + 2 gi 3 [/xiXiPll + fxsx^Psi + ‘^fxsXiPllPSl], 

F^2 = 2gii [f^,^^PllPl2 + fx2X2P2lP22 + fx3X3P3lP32 + fx^xi {PllP32 + Pl2P3l) + fxiX2 {Pl2P21 + 

P 22 PI 1 )] + 2gi2[/Jia:iPllPl2 + fx 2 X 2 P 2 lP 22 + fxiX 2 iPl 2 P 21 + P 22 Pn)] + ‘2qi3[fxixiPim2 
+fx3X3P3lP32 + fx 3 Xi {PnP32 + PuPSl)], 

F 22 = 2gil + fl2X2P22 + fx3X3P32 + ^fx3X3Pl2P32 + “2 fl^^^Pl2P22] + 2^12 [/Jia;iPl2 

+fx 2 X 2 P 22 + ‘ 2 fxiX 2 P 32 P 22 ] + 2^13 [/|ixiPl2 + + fx 3 XiPl 2 P 32 )], 

Fhl = 6gii6ii [f^^^,PllPl3 + fx2X2P2lP23 + fx3X3P3lP33 + fx3Xi iPllP33 + Pl3P3l) + fx3X2 iPllP23 
+P 21 P 13 )] + Gqi2bn[fx3XiPnPl3 + fx2X2P2lP23 + fx3X3P3lP33 + fxiX2iP33P23 + Pl3P2l)] 
+Qqi3bn[fxixiPim3 + fxsxsPsms + fx3Xi{PnP33 + Pl3P3l)], 

^122 = 2^11 {2pi2Pl3bl2 + PIIP 13 & 22 ) + fx2X2 {2P22P33bl2 + P21P23&22) + (2p32P33^12 

+P3lP33b22) + fx3Xii2pi3P32bl2 +PllP33^22 +Pl3P31^22 + 2^12^33^12) + (2^12^23 

x6i 2 +P13P21&22 +PllP23^22 + 2p22Pl3^12)] + 2gi2 (2pi2Pl3fel2 +PllPl3^22) + fx2X2 

x(2p22P23&12 +^ 21 ^ 23 ^ 22 ) + fxiX2i‘^P32P23bl2 + PllP23^22 + 2^13^22^12 +P 13 P 2 I& 22 )] 

+ 2 gi 3 [/i)ia;i( 2 pi 3 Pl 2 fel 2 +PllPl 3 ^ 22 ) + /^axa ( 2 P 32 P 33 bl 2 + P 31 P 33 & 22 ) + ( 2 ^ 32 ^ 13^12 

+P 3 lPl 3 fe 22 + 2 pi 2 Pl 3^12 +P33PII&22)], 
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F?i 


Fh 


F ^2 = 


Fh 2 


Tp 2 

-^222 


^Q2l[fxiXiPll 1x2X2^21 fx^x^Psi “^fx^xiP^^P^^ ‘^fxiX2P^^P‘^^^ ‘^Q22[fxiXiPll /; 

xpil + 2/Ji^2PllP2l] + 2g23[/xixiPll + fxsx^pll + ‘^fxsXiPllPSl], 

^q2l[fxixiPim2 + fx2X2P2lP22 + fx3X3P3lP32 + fx3XiiPllP32 + PuPSl) + /iiX2 ^12^21 + 
P22Pn)] + ‘^q22[fxixiPnPl2 + fx2X2P2lP22 + /xi^z ^12^21 + P22Pn)] + 2 ^ 23 [/xia;iPllPl2 
+fx 3 X 3 P 3 lP 32 + fx3Xi {PnP 32 + Pl 2 P 3 l)], 

‘^Q 2 l[fX1X1P12 F fx 2 X 2 P 22 fx 3 X 3 P 32 F ‘^fx 3 XiP 32 P 32 “1“ 2/2,^^3,2^12^22] “1“ 2(^22[/a;ia;iPl2 

+/x 2X2P22 + 2/Jia;2^^12P22] + 2^23 [/xixiPl2 + + /x 3 XiPl 2 P 32 )], 

292l[/iixi(2PllPl3&12 +Pl2Pl3^1l) + fx2X2i^P2lP23h2 +P22P23&ll) + /^xa (2?»3lP33^12 
+P 32 P 33 fell) + /i 3 a;i( 2 pilP 33&12 +Pl 3 P 32 bll +P 12 P 33&11 + 2 ^ 13 ^ 31 ^ 12 ) + /iix2 ( 2 Pl 3 P 21 
x6i2 +P 12P23&II + 2pilP23fel2 +P22Pl3^1l)] + 2(/22 [/Jixi (2pilPl3fel2 +Pl2Pl3^1l) + fx2X 

xi2p2lP23bl2 F P22P23bn) F /JiX2(2PllP23^12 +P12P23&11 + 2pi3P21&12 + Pl3P22fell)] 

+2g23[/xixi(2pilPl3&12 +?»12Pl3^1l) + /xa^a (2p3lP33^12 +P32P33^1l) + /xgxi (2p3lPl3^12 
+P 33 Pl 2 ftll + 2P11P33612 +P 32 Pl 3 bll)], 

6921 ^ 22 [/iixiPl2Pl3 + fx2X2P22P23 + fx3X3P32P33 + fx3XiiPl2P33 F P 13 P 32 ) + fx3X2iPl2P23 
+P 13 P 22 )] + 6Q'22&22[/sia:iPl2Pl3 + fx2X2P'i2P23 F /xiX2 (^*12^23 +P 13 P 22 )] 

+6^23^22 [/x, XIP 12 P 13 + fLx^P32P33 + fX'liXl (Pl2P33 F P13P32)]- 


•2 

X2X2 


6.2 Global stability of the system (2.1) around i?* 

Theorem 6 . 1 . The interior equilibrium is globally asymptotically stable if the condition 
(i) aia2bib2rk > {xu F k){w + k){aibiC2 + a2^2Ci) is satisfied. 


Proof. Let 


L(xi,a; 2 ,a; 3 ) = Li(xi, 2:2, a^s) + L 2 (xi,X 2 ,X 3 ) + L 3 (xi,X 2 ,X 3 ) (6.1) 

be a positive Lyapunov function, where 

Li = si(xi - Xu - Xu ln(—)), L2 = S2{x2 - X2* - X2* ln(—)), 

3 ^ 1 * ^ 2 * 

L3 = 53(2:3 - 2:3* - 3:3* ln( —)); 

X3* 

si, S2 and S3 being positive real constants. 

This function is well-defined and continuous in Int(i 2 +^). It can be easily verified that the 
function L(xi,X2,a:3) is zero at the equilibrium point E^, and is positive for all other positive 
values of (xi,X2,X3), and thus E^ is the global minimum of L(xi,X2,X3). 

Since the solutions of the system are bounded and ultimately enter the set 11 = {(xi, X2,2:3); xi > 
0 , X2 > 0 , X3 >0 ■ xi + ^ + ^ < M+e, V e > 0 }, we restrict our study in 11 . The time derivative 


13 


of L along with the solutions of the system (2.1) gives (cf. Sarwardi et al. [25], [29]) 


dt 


X 

+ 

X 


-Si 


rk 


-(xi* + k){xi + k) 


C 1 X 2 * 

(ai + xi + 6 iX2)(oi + Xu + biX2*) 


^2^3* / \2 

(02 + Xi + 62 X 3 )(a 2 + Xi* + 62 X 3 *) J 

/ _ s 2 _ [_ € 262(02 + ^ 2 ) _ 

X2 X2* '^3 

€1(526161X2* - Si(ai + Xi*))(xi - X1*)(X2 - X2*) 


(oi + Xi + biX2){ai + Xu + 61X2*) 
(xi - X1*)(X3 - X3*). 


+ 


■_ €i6i(ai + 61) _■ 

-(oi + Xi + 6iX2)(ai + Xi* + 61X2*)- 

(X 3 - X 3 *)^ 

62(536262X3* - Sl(o2 + Xl*)) 

(02 + Xl + 62X2)(a2 + Xl* + 62X2*) 

( 6 . 2 ) 


Letting 51 = 1, 52 = and 53 = we have 


dL 

dt 


< - 


rk 


€ 1 X 2 * 


-(xi* + A:)(xi + k) (ai + xi + 6 iX 2 )(ai + xi* + 61 X 2 *) 


62 X 3 * 


< - 


(02 + Xl + 62X3) (02 + Xl* + 62X3*)- 
rk Cl €2 ■ 


.{xu + k){w + k) aibi 0262 - 
< 0, by condition (i), 


(xi - Xl*)^ 
(x — x*)^ 


(6.3) 


along all the trajectories in the positive octant except (xi*,X2*,X3*). Also ^ = 0 when 
(xi,X2,X3) = (xi*,X2*,X3*). The proof follows from (6.1) and Lyapunov-Lasalles invariance 
principle (cf. Hale [30]). 


Table 1: Schematic representation of our analytical findings: LAS = Locally asymptotically 
stable, GAS = Globally asymptotically stable, HB = Hopf bifurcation, SHB = Subcritical Hopf 
bifurcation. 


Equilibria 

Feasibility conditions/ 
parametric restrictions 

Stability conditions/ 
parametric restrictions 

Nature 

Eo 

No Gonditions 

No Conditions 

Unstable 

El 

6161 > 61, Xl, > 

6iXi,+X2, 

LAS 

E2 

6262 > 62, XI2 > 

62XI2 + X 32 > ^ " 

LAS 


Xl* 

A: < min{ai + 61x2*, 02 + 62x3*} 

LAS 

E* 


(i) r > 61 + 62, (n) Xl, > 

, (in) Xl, > 

Persistence 

E* 


Stated in the Proposition 3.3 

Boundedness 

E* 


n > 0 (cf. equation (6.15)) 

SHB 

E* 


aia2bib2rk > {xu+k){'w+k)x 
(016162 + 026261) 

GAS 
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Table 2: The set of system parameter (including the critical parameter rc) values and their 
corresponding Figures with description. 


No. 

Fixed Parameters 

r 

Figures 

Description 

1 

r = 1.37 >rc = 1.320961640, k = 
200 , ai = 100,02 = 100,61 = 
0 . 5,62 = 0.5; Cl = 1.8, C 2 = 
1 . 8 ,5i = 0.82,52 = 0.62, Cl = 
0.8143,62 = 0.6250.. 

1.37 

Figs. 1: 
(a)-(b) 

2D view of Hopf 
bifurcation 

2 


1.320961640 

Figs. 2 

Limit cycle 

3 


r G [0.8, 2.0] 

Fig. 3 

Hopf bifurcation 

(growth rate r vs. 
population volumes) 

4 


1.4700000000 

Fig. 4 

2D view of local 
stability 

5 


1.4700000000 

Fig. 5 

3D view of local 
stability 

6 


1.4700000000 

Fig. 6 

Global stability 

7 

r = 1.37, A: = 200, ai = 100,02 = 
100,61 = 0.5; 62 = 0.5, Cl = 
1 . 8 , c 2 = 1 . 8 ,5i = 0 . 82,52 = 
0.62,61 = 0.8143,62 = 0.6250 

rc = 1.320961640, 
H = 1.0424314050 

Figs. 7 

Subcritical Hopf 
bifurcation 


7 Numerical simulation 

For the purpose of making qualitative analysis of the present study, numerical simulations have 
been carried out by making use of MATLAB-R2010a and Maple-12. The analytical findings of 
the present study are summarized and represented schematically in Table 1. These results are 
all verified by means of numerical illustrations of which some chosen ones are exhibited in the 
figures. Here, we have given some numerical simulations on the study of stability and bifurcation 
of the proposed system (2.1) around the interior equilibrium E^. We took a set of admissible 
parameter values: r = 1.7, A: = 200, oi = 02 = 100, 61 = 62 = 0.5, ci = C 2 = 1.8,= 0 . 82,52 = 
0.62, ei = 0.8143 ,62 = 0.6250. For this set of parameter values, it is found that the system 
possessed an unique interior equilibrium point = (169.1663564,55.36073780,62.98120968). 
The system parameter r is the growth rate of the prey population which plays a crucial role 
in regulating the dynamical behaviour of the proposed system. For this reason, we take this 
parameter as an influential parameter and try to determine the possible outcomes by varying 
this parameter within its feasible range. The interior equilibrium E^ is stable for the values 
of r > Tc = 1.320961640 (cf. Figure: 4-5 for local stability and Figure: 6 for global stability). 
The system (2.1) experiences Hopf bifurcation when the parameter r crosses the critical value 
Tc from left to right, i.e., when r = Tc, all the species coexist in the form of periodic oscillation. 
Following the steps discussed in Subsection 6.1, we have found the value of H = 1.042431405 > 0, 
which indicates that the the obtained Hopf bifurcation is subcritical bifurcation (cf. Figure 7). 

It is observed that, if the interference coefficient hi (interference effect due to the presence of 
second predator on the first predator) increases it stabilize the system for hi = 0.7 while it 
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is unstable at bi = 0.6 and when the parameter 6 i exceeds its value 20 , the first predator 
population is died out from the system, i.e., the system breakdown. Similarly, the interference 
effect due to the presence of first predator on the second predator, parameterized by 62 plays 
an important role to stabilize the system. If the interference coefficient 62 increases it stabilizes 
the system for 62 = 0.6, while it is unstable at 62 = 0.5. It also regulates the existence of 
second predator in the system. As the parameter 62 exceeds its value 20.9, the second predator 
population is died out from the system. 

Analogously, if the parameter di, denoting the death rate of first predator increases then the 
volume of the fist predator decreases as well as second predator population increases and if di 
decreases, the first predator population increases and second predator population decreases. If 
the death rate di is gradually increased to a certain level the first predator population goes 
into extinction. Similar result is observed for the case of the second predator’s death rate. 
The above observations ensure that the model under consideration is consistent with biological 
phenomenon (Figures are not reported here). 



Figure 1: 2D view of Hopf bifurcation around the interior equilibrium A* of the system (2.1) with 
parameter values: r = 1.37 > Vc = 1.320961640, fc = 200, ai = 100,02 = 100, 61 = 0.5, &2 = 0.5; ci = 
1.8, C 2 = 1.8, 5i = 0.82 ,52 = 0.62, ei = 0.8143 ,62 = 0.6250. 


8 Concluding remarks 

The problem describes by the system (2.1) is well posed that xi, X 2 and X 3 axes are invariant 
under the flow of the system. So far our knowledge goes this is the first attempt to study an 
ecological system with semilinear/bilinear growth of the prey population. Generally, researcher 
only studied biological model systems with logistic/linear growth of prey population. Here is the 
novelty of our study. One of the important observations is that the prey population becomes 
unbounded in absence of its admissible predator in long run of time. But in the presence 
of predator species the prey population can be made bounded under suitable combination of 
system parameters and as a consequence it is shown that the total environmental population 
under consideration is bounded above (cf. Subsection 3.3). Therefore, any solution starting in 
the interior of the first octant never leaves it. This mathematical fact is consistent with the 
biological interpretation of the system. Due to the inclusion of semilinear/bilinear growth of the 
prey population, the axial equilibrium point is driven away by the system, which is rarely found 
in the modern research work on Mathematical biology. Thus, the prey population alone can not 
survive in stable condition without their admissible predator populations. It is found that only 
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Figure 2; Limit cycle behaviour of the dynamical system at -E* with the same parameter values 
used for Figure 1. 


the mutual interference between the predators, which are parameterized by bi and 62 can alone 
able to stabilize the prey-predator interactions even when a semilinear/bilinear intrinsic growth 
rate of prey population is considered in the proposed mathematical model. Whereas these 
parameters have much contribution in stabilizing prey-predator interactions when only linear 
intrinsic growth rate is considered in some mathematical models (cf. Dimitrov and Kojouharov 
[31]). It is observed in the study of this model system that there exist a balance between the 
predator’s need for food and its saturation level and in this case is likely to be expect a periodic 
behaviour in long run. This behaviour is neutrally stable but relatively unstable. A small change 
in the parameters (caused by environmental changes for instances) forces the system to stabilize 
around the interior equilibrium or to oscillate indefinitely around interior equilibrium (by going 
away from it, which causes collapse of the system or breaks the coexistence of the population). 
Representative numerical simulations of this case are shown in Figures: 1-3, which support our 
analytical findings (cf. Theorems 5.2 and 6.1). We have also established the sufficient conditions 
for the global stability of the coexistence equilibrium (cf. Figures: 5-6). 
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Ph.D. supervisor Prof. Prashanta Kumar Mandal, Department of Mathematics, Visva-Bharati 
(a Central University) for his generous help while preparing this manuscript. 
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asymptotic stable. Here the same set of parameter values is used for Figure 4 except r = 1.29. 
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